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Abstract 

This  document  represents  the  final  report  on  the  various  scientific  activities  and  ac¬ 
complishments  relating  to  Grant  No.  FA9550-08-1-0151  over  its  period  of  performance, 
April  1,  2008  -  July  31,  2011. 

The  accomplishments  may  be  subdivided  according  to  the  project’s  theoretical, 
experimental,  and  post-processing/computational  objectives.  Among  the  theoretical 
accomplishments,  we  list 

•  a  new  model  for  a  sparse  representation  of  man-made  space  objects  and  its  use, 
via  a  new  spectral-correlation  approach,  to  segment  their  material  components; 

•  the  use  of  2D  segment-boundary  data  for  multiple  poses  of  a  man-made  space 
object  to  recover  its  3D  shape; 

•  new  fundamental  results  involving  statistical  information  and  Bayesian  error 
bounds  for  characterizing  the  performance  of  reconstruction  algorithms,  and  feature- 
extraction  and  estimation  fidelity;  and 

•  the  use  of  Fisher  information  and  the  associated  Cramer- Rao  lower  bound  on 
estimator  variance  to  characterize  the  value  of  a  prior  knowledge  of  object  support 
for  spatial-bandwidth  extension  beyond  diffraction-limited  observations. 
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The  experimental/instrumentation- related  accomplishments  include 

•  new  CASSI  instrumentation  for  supplying  multiple  frames  of  high-resolution, 
broad-band  hyperspectral  image  (HSI)  data  in  compressive  form; 

•  supplying  CS  data  from  CASSI  to  the  WFU-based  researchers  to  test  their  algo¬ 
rithms  for  efficiency,  accuracy,  and  robustness  of  HSI  reconstructions  from  single¬ 
frame  and  multi-frame  data;  and 

•  developing  a  push-broom  variant  of  CASSI  for  comparing  ground  truth  HSI  data 
with  reconstructions  from  compressively  sensed  (CS)  CASSI  data. 

Finally,  the  computational/mathematical  accomplishments  include 

•  developing  and  testing  a  number  of  new  joint  segmentation  and  spectral  recovery 
algorithms  for  uncompressed  HSI  data; 

•  developing  an  alternating  LS  algorithm  for  segmenting  and  reconstructing  the  full 
HSI  datacube  from  compressed  CASSI  data  and  testing  it  on  both  simulated  and 
real  experimental  data;  and 

•  a  sparse  nonnegative  matrix  under-approximation  approach  that  uses  a  sparsity 
constraint  to  process  HSI  data. 
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1  Overall  Project  Objectives 

To  place  the  project  accomplishments  in  proper  perspective,  we  first  present  its  originally 
proposed  objectives.  Specifically,  we  aimed 

1.  to  develop  fundamental  mathematical  models  and  system  design  methodologies  sup¬ 
porting  the  identification  of  space  objects  from  multi-modal  data; 

2.  to  simulate  and  design  laboratory  experiments  to  demonstrate  novel  compressive  spec¬ 
tral  coherence  imaging  approaches  for  SOI; 

3.  to  develop  an  information-theoretic  toolbox  of  concepts  and  metrics  to  analyze  system 
performance;  and 

4.  to  offer  numerical  approaches  that  yield  the  most  computationally  efficient  algorithms 
for  post-processing  the  data  generated  by  those  methods. 

The  project  was  subdivided  according  to  its  theoretical,  experimental,  and  mathemati¬ 
cal/computational  thrusts,  with  the  three  collaborating  institutions,  UNM,  Duke,  and  WFU, 
being  the  respective  leaders  of  these  three  thrusts.  We  begin  with  a  summary  of  UNM’s  work 
in  meeting  the  theoretical  objectives  of  this  effort,  followed  by  Duke  U.’s  work  leading  to 
an  experimental  realization  of  compressive  spectral  sensing  in  the  presence  of  atmospheric 
turbulence,  and  end  with  WFU’s  contributions  to  algorithm  development  for  various  com¬ 
pressive  sensing  reconstruction  strategies.  The  work  in  these  three  research  thrust  areas  has 
resulted  in  a  series  of  important  publications  and  invited/contributed  presentations  at  pro¬ 
fessional  meetings  in  computational  imaging  and  data  analysis.  These  are  listed  in  Sections 
3-5  of  this  report. 

2  Summary  of  Overall  Project  Accomplishments 

The  four  originally  proposed  objectives  of  the  effort  were  largely  achieved  over  the  project 
duration.  In  this  report  we  provide  only  a  summary  description  of  the  project  accomplish¬ 
ments,  leaving  the  many  technical  details  of  a  fuller  description  to  the  project-supported 
publications  that  are  cited  at  appropriate  places  below  and  are  publicly  accessible.  Our 
technical  accomplishments  may  be  summarized  as  follows: 

2.1  Sparse  Hyperspectral  Object  Model 

We  proposed  and  exploited  a  unique  model  of  sparsity  of  satellite  objects.  Such  objects  may 
be  regarded  quite  accurately  as  being  composed  of  homogeneous  material  segments  with 
uniform  spectral  and  brightness  characteristics  except  possibly  near  the  segment  boundaries. 
This  implies  a  relatively  small  number  of  characteristic  parameters  when  compared  with  the 
vast  number  of  spectral-spatial  voxels  in  the  hyperspectral  datacubes  of  such  objects.  This 
observation  was  exploited  to  recover  the  values  of  these  parameters  and  thus  reconstruct  the 
whole  hyperspectral  datacube  from  both  compressive  and  non-compressive  noisy  spectral 
data.  A  number  of  papers  published  by  the  project  team  chronicle  this  development  over 
the  project  duration. 
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2.2  A  Spectral  Correlation  Approach  for  Segmentation  and  Spec¬ 
tral  Recovery 

Hyperspectral  datasets  typically  contain  far  more  information  than  needed  for  a  typical 
remote  sensing  task  of  material  identification,  segmentation,  and  composition.  A  typical 
solar-illuminated  man-made  satellite  consists  of  a  relatively  small  number  of  geometrical 
primitives  and  materials  with  unique  spectral  reflectance  signatures.  This  implies  a  highly 
correlated  hyperspectral  data  cube,  which  may  be  idealized  as  a  low-rank  sum  of  terms, 
each  expressible  in  a  spatial-spectral  factorized  form.  The  high  degree  of  spectral  correlation 
within  each  spatial  segment  of  the  object  admits  a  simple  moving-averages  (MA)  approach 
for  identifying  the  material  components  and  their  spatial  boundaries,  as  we  have  recently 
demonstrated  (Li,  Ng,  Plemmons,  Prasad  and  Zhang,  2010).  This  approach  works  equally 
well  with  a  less  idealized  form  of  the  data  cube  that  allows  for  arbitrary,  non-uniform  values 
within  the  support  of  each  material  component. 

In  the  MA  approach,  we  cross-correlate  the  spectral  traces  at  different  spatial  pixels  in 
the  dataset  with  the  spectral  traces  of  different  materials,  either  picked  from  a  sufficiently 
exhaustive  database  of  spectra  or  extracted  from  the  HSI  dataset  itself.  To  improve  the  SNR 
of  the  normalized  cross-correlation,  we  average  it  over  a  small  cell,  say  2x2  or  4x4  square 
pixels  in  size,  that  we  move  over  the  full  spatial-pixel  array  in  the  HSI  dataset.  As  the 
averaging  cell  moves,  the  normalized  correlation  varies  between  0  and  1.  The  value  1  signals 
a  perfect  match  between  the  test  trace  and  the  spectral  trace  of  pixels  inside  the  moving  cell, 
while  a  much  smaller  value  indicates  poorly  matched  or  mismatched  spectra.  By  repeating 
this  calculation  with  different  test  traces,  we  can  then  segment  the  various  components  and 
also  identify  their  spectra. 

To  demonstrate  this  approach,  we  used  the  hyperspectral  image  (HSI)  data  cube  of  the 
Hubble  space  telescope  (HST)  that  was  simulated  in  Zhang,  Wang,  Plemmons,  and  Pauca 
(2010)  for  demonstrating  low-rank  tensor  methods  on  HSI  data.  Figure  1  shows  a  noisy  HST 
spatial  brightness  distribution  (PSNR=20)  in  the  mid-spectral  band.  The  effectiveness  of 
the  moving-averages  approach  in  segmenting  the  eight  different  material  components  of  the 
HST  from  its  HSI  data  cube  is  easily  seen  in  the  eight  panels  of  Fig.  2. 

2.3  Surface  Shape  Recovery  and  Pose  Estimation 

The  unique  relationship  among  the  two-dimensional  (2D)  projections  of  a  curve  fixed  on  a 
rigidly  rotating  three-dimensional  (3D)  object,  with  one  projection  per  pose,  implies  that 
those  projections  contain  precise  information  about  the  object  shape  and  size  parameters  as 
well  as  the  Euler  angles  describing  the  poses  assumed  by  the  object  as  it  rotates.  We  can 
exploit  this  observation  for  a  rotating  3D  space  object  by  following  the  kinematical  evolution 
of  its  segment  boundaries  whose  2D  image-plane  projections  can  be  determined,  e.g.,  by  the 
moving-averages  approach.  If  each  material  component  of  the  object  has  a  specific  but 
simple  geometry,  we  can  idealize  the  components  in  terms  of  simple  geometrical  primitive 
shapes,  such  as  spheres,  ellipsoids,  cylinders,  cones,  parallelepiped,  triangles,  tetrahedra, 
etc,  drawn  from  a  database  of  shapes.  Given  the  2D  image-plane  projections  of  a  curve 
and  a  certain  primitive  shape,  we  construct  a  fit-to-data  L2  functional  to  optimize  for  the 
shape  parameters,  Euler  angles  for  each  pose,  the  3D  body  coordinates  of  points  on  the 
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Figure  1:  The  spatial  brightness  distribution  of  the  simulated  HST  in  the  middle  spectral 
channel 


Figure  2:  Map  of  the  spectral  cross-correlation  function  for  the  eight  different  material  traces 
used  in  the  construction  of  the  HST  datacube. 


curve,  and  the  rotation  matrix  for  each  pose.  For  this  highly  nonlinear  problem,  we  use  a 
block  coordinate  descent  algorithm  to  alternately  optimize  for  each  set  of  parameters  in  an 
iterative  way.  And  by  searching  exhaustively  through  a  database  of  primitives,  we  can  infer 
the  specific  primitive  shape  and  its  size  parameters  that  correspond  most  closely  to  the  time 
series  of  observed  2D  projections  of  a  specific  boundary. 

In  our  initial  studies  (Prasad  and  Zhang,  2010),  through  forward  models  we  simulated  the 
2D  projections  of  curves  on  an  ellipsoid  and  a  cylinder  at  two  different  poses,  from  which  we 
sought  to  recover  the  original  parameters  and  body  coordinates  used  in  the  simulation.  The 
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Figure  3:  Two  curves  on  an  ellipsoid  of  certain  orientation,  shown  in  red,  fitted  to  an  ellipsoid 
via  an  iterative  fitting  algorithm 


results  have  been  very  promising.  We  have  shown,  in  particular,  that  with  just  a  few  poses 
of  the  rotating  surface  and  its  segment  boundaries,  the  fit  to  wrong  shapes  across  segment 
boundaries  can  be  heavily  discriminated  against.  Even  with  just  two  poses,  we  could  infer 
the  correct  shape  from  similar  looking  projections  while  discriminating  against  other  shapes, 
as  we  see  in  Figs.  3  and  4.  Here  the  two  poses  of  a  curve  on  an  ellipsoid  are  shown  in 
red.  In  Fig.  3  we  fitted  the  two  poses  to  a  correct  ellipsoidal  surface  shape,  and  determined 
the  parameters  of  the  ellipsoid  iteratively.  One  can  see  that  after  about  10  iterations  of  the 
fitting  alorithm,  the  fit,  shown  by  blue  curves,  is  quite  good.  By  contrast,  when  the  curves 
are  fitted  to  a  cylinder,  which  is  a  wrong  shape,  the  fit  does  not  improve  at  all  over  the  first 
20  iterations  and  beyond,  as  shown  in  Fig.  4.  Estimated  shape  parameters  and  Euler  angles, 
not  given  here,  are  also  quite  close  to  the  original  when  the  correct  shape  was  tried. 

This  approach  can  be  exploited  to  devise  a  recursive  algorithm  for  determining  the  surface 
geometry  of  the  entire  object  as  it  rotates  to  reveal  its  complete  surface  over  a  certain 
observation  period. 

Surface  shapes  can  be  further  constrained  by  well  resolved  polarimetric  observations  of 
the  reflected  sunlight.  Under  recently  obligated  AFOSR  funding,  we  are  now  developing, 
testing,  and  extending  our  technique  further  to  include  spectro-polarimetric  image  data  in 
order  to  achieve  a  more  comprehensive  characterization  of  the  surface  properties  via  purely 
optical  observations.  We  are  also  exploring  the  use  of  superquadric  parametric  surfaces  to 
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Figure  4:  Two  curves  on  an  ellipsoid  of  certain  orientation,  shown  in  red,  fitted  to  an  ellipsoid 
via  an  iterative  fitting  algorithm 
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simulate  more  general  surfaces  than  our  simple  primitive  shapes  have  to  offer. 

2.4  Fundamental  Relations  between  Statistical  Information  and 
Bayesian  Error  Bounds 

A  number  of  mathematical  relations  connecting  statistical  information  (SI)  with  Bayesian 
error  analysis  were  developed  and  analyzed  under  the  grant  effort.  In  particular,  we  have 
studied  just  how  the  minimum  mean-squared  error  (MMSE)  of  Bayesian  estimation  and 
minimum  probability  of  error  (MPE)  in  Bayesian  M-ary  hypothesis  testing  bound  the  equiv¬ 
ocation  entropy,  which  represents  the  loss  of  information  in  a  communication  channel  such 
as  an  imaging  system.  These  quantities  were  calculated  for  the  problem  of  localizing  an  un¬ 
resolved  source  with  sub- diffractive  error  in  the  presence  of  image  blur  and  detection  noise. 
This  calculation  is  evidently  quite  relevant  for  locating  faint  space  objects,  such  as  space  de¬ 
bris  and  small  GEO  satellites  that  are  unresolved  in  the  space  imaging  setting.  The  details 
of  this  work  may  be  found  in  Prasad  (2011)  (submitted  to  JOSA  A),  Prasad  and  Narravula, 
2011,  and  some  unpublished  work  (Narravula  and  Prasad);  here  we  present  only  important 
results. 


Relations  Involving  MPE,  Statistical  Equivocation,  and  Mutual  Information  An 

important  result  in  this  context  is  the  following  that  relates  the  MPE  for  M-ary  detection  to 
a  generalized  version  of  the  equivocation  entropy  (GEE)  of  statistical  information  theory: 


p(min)  _  1  _  E 


exp 


(°°  =  *) 


(a 


where  the  GEE  of  nth  order,  Hn{V.\X  =  x),  is  defined  as 


M 


Hn(H\X  =  x)  =  Pn(m\x)  In  Pn(m\x), 


(2) 


m= 1 


in  terms  of  the  following  n-indexed  PDF  for  a  specified  data  value  x: 

n  /  in  Pn(m\x) 

Pn(m\x  =  ,  m  =  1, . . . ,  M. 

EmP  (mF) 


(3) 


From  this  exact  result,  there  follow  a  number  of  well  known  bounds  on  the  MPE,  namely 
the  Santhi-Vardy  and  the  Renyi  bounds  involving  the  ordinary  equivocation  entropy  related 
more  directly  to  the  MI. 

The  MI  can  be  shown  to  be  bounded  below  by  a  logarithmic  function  of  MMSE,  as  shown 
by  Seidlcr  nearly  forty  years  ago.  However,  there  is  also  an  upper  bound  related  to  the  MPE, 
given  data  X,  that  we  have  derived  recently,  namely 


/(©;  X)  <H(Q)+EX 


ma x{p(9m\X)}  logrnax{p(6»m|X)}  , 


(4) 


where  H(Q)  is  the  signal  entropy.  The  tightness  of  the  upper  bound  (4)  is  directly  determined 
by  the  smallness  of  the  average  amount  of  spill-over  of  any  posterior  probability  mass  function 


Minimum  Error  Probability  for  M-fold  OSR  with  Poisson  Data  MMSE  for  M-fold  OSR  with  Poisson  Data 


L/eo  L/eo 


(a)  (b) 

Figure  5:  (a)  Minimum  MPE  and  (b)  MMSE  for  resolving  the  source  position  vs.  image 
pixel  width  (in  units  of  diffraction-limited  pixel  width)  for  M  —  4  for  three  different  source 
strengths,  K  =  10,  103,  and  105,  in  the  Poisson  shot-noise  model. 


(PMF),  p(0m\X),  into  neighboring  decision  regions,  namely  Imrn/ ,  m'  ^  m,  that  we  can  also 
use  to  relate  the  MPE  to  MMSE  (see  next  paragraph).  For  this  reason,  this  bound  is  expected 
to  be  quite  tight  for  highly  sensitive  Bayesian  detectors,  A  number  of  other  new  bounds  on 
the  MI  in  terms  of  the  MPE  have  also  been  obtained,  but  because  of  space  considerations 
we  are  not  presenting  them  here. 

Relation  between  MPE  and  MMSE  and  Application  to  Source  Location  For  a 

highly  sensitive  detector,  i.e.,  when  the  SNR  is  high,  the  MMSE  and  MPE  are  both  related 
directly  to  the  probability  for  incorrect  inference,  namely  that  one  of  the  incorrect  M  —  1 
parameter  values  gave  rise  to  the  data.  If  Imm>  denotes  the  probability  of  inferring  the 
parameter  value  9m>,  given  the  parameter  value  6m,  m  ^  rn' ,  then  the  MMSE  and  MPE  are 
both  related  to  this  quantity  as 

MMSE  --  EE  {dm  drn/')  {Immf  dr  ^ m'm )  dr  0(6  ) 

m  m' 

p,„„n,  =  1^12  +  W  •  (*>) 

m  m'y^m 

where  terms  of  smaller  order  can  be  neglected  in  the  first  expression  in  the  limit  of  large 
SNR.  This  close  connection  between  the  MPE  and  MMSE  is  seen  in  Fig.  5  for  a  single-pixel 
point-source  localizer  for  photon-limited  imaging. 

We  have  derived  a  tight  upper  bound  on  the  MMSE  (Narravula  and  Prasad,  unpublished), 
which  makes  it  possible  to  calculate  the  MMSE  approximately  in  the  most  general  situation 
thus  rendering  it  numerically  highly  amenable.  This  bound, 

MMSE  <  E[59(k)  -  Sd}2,  (6) 
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Figure  6:  Plot  of  the  source-strength  threshold,  Kmin,  vs.  the  degree  of  localization,  M,  for 
the  one-pixel  imager  (square  symbols)  and  the  three-pixel  imager  (+  symbols). 


where  the  prefix  5  denotes  linear  deviations,  8  A  =  A  —  E(A),  and  6  is  any  estimator  of  6,  can 
be  made  tight  by  a  variational  approach  in  which  one  minimizes  the  RHS  of  this  inequality 
over  a  parametric  family  of  test  estimators.  Use  of  this  bound  enables  us  to  derive  the 
photon-number  thresholds  needed  to  achieve  point-source  localization  with  sub-diffractive 
error,  a  type  of  optical  super-resolution.  The  minimum  photon  number  needed  for  M-fold 
reduction  in  localization  error  relative  to  diffractive  localization  has  an  expected  M2  form 
when  several  pixels  worth  of  image  data  are  exploited  for  this  purpose,  but  a  single  pixel 
image  requires  a  far  steeper  photon  threshold  of  form  M2  77,  as  displayed  in  Fig.  6  below. 

2.5  Statistical  Information  in  Reconstruction  Algorithms 

Another  application  of  our  SI  based  performance  metrics  has  involved  the  study  of  informa¬ 
tion  recovered  by  an  algorithm  that  exploits  image  data  with  or  without  any  post-processing 
constraints  such  as  the  knowledge  of  object  support.  We  have  generated  plots  of  mutual  in¬ 
formation  (MI),  namely  the  information  in  the  ensemble  of  object  scenes  that  is  successfully 
transmitted  by  an  imaging  system  to  its  reconstruction  algorithm,  with  respect  to  varying 
amounts  of  detection  noise  and  object  support  sizes  to  quantify  information  theoretic  perfor¬ 
mance  upper  bounds.  The  imaging  system  of  interest  here  has  been  one  based  on  compressed 
sensing  (CS)  that  performs  projective  measurements  of  masked  image  scenes  involving  two 
or  more  homogeneous  surface  components  with  different  brightness  values  and  then  uses  an 
algorithm  to  extract  the  position,  size,  and  material  signature  (via  pixel  brightness  differ¬ 
ences  between  neighbors)  of  the  components.  We  have  showed  that  the  MI  obtained  at  the 
post-processing  (algorithm)  stage  is  lower  than  at  the  data  stage,  since  any  practical  algo¬ 
rithmic  reconstruction  has  finite  errors  of  its  own.  However,  the  use  of  a  priori  knowledge, 
such  as  constraining  the  recovery  of  the  center  position  of  a  component  to  a  pre-specihed 
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MI(X,X)  vs.  noise  level 


Figure  7:  Plot  of  the  mutual  information  (MI)  vs.  detection  noise  for  the  estimator  algorithm, 
as  compared  to  the  MI  in  the  CS  data  themselves. 


rectangular  support  region  about  the  true  center,  can  be  shown  to  improve  the  MI  beyond 
that  contained  in  the  data  alone  if  the  support  dimensions  are  chosen  to  be  sufficiently  small. 

The  following  lower  bound  on  the  MI  codifies  this  relationship  between  the  MI,  I(X ;  X), 
and  the  variance  of  any  estimator  A": 

/(A;  X)  >  h( X)  -  h(X\X)  <  ^E£.  In  [27reEX|* (x  -  A)2]  ,  (7) 

where  h( A)  is  the  differential  entropy  of  the  object,  which  is  the  maximum  information  that 
a  perfect  noiseless  measurement  can  achieve.  The  fact  that  the  estimator-recovered  MI  is 
smaller,  in  general,  than  the  data-contained  MI,  /(A;  Y),  is  borne  out  in  Fig.  7,  where  we 
plot  the  MI  vs.  detection  noise  for  a  readout-noise  dominated  compressive  sensing  system 
that  extracts  the  center  position  of  a  square  of  given  size  inside  another,  larger  square 
background.  Note  that  by  sharpening  the  prior  constraint,  in  this  case  the  support  of  the 
region  of  uncertainty  of  the  center  position  of  the  inner  square,  we  can  more  than  overcome 
the  loss  of  information  caused  by  the  detection  noise  when  that  is  large. 

2.6  Fisher-Information  Analysis  of  Support- Assisted  Optical  and 
Digital  Superresolution 

Early  in  the  project,  we  expolited  Fisher  information  (FI)  to  characterize  the  extent  of 
spatial-frequency  extrapolation  beyond  the  diffraction-limited  optical  bandwidth  when  the 
support  of  the  object  is  known  a  priori.  This  support-assisted  optical  superresolution  (OSR) 
phenomenon  was  analyzed  for  both  ID  and  2D  images  by  means  of  a  Shannon- Nyquist  type 
of  sampling  expansion  in  the  spatial-frequency  domain  for  a  spatially  limited  object. 
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A  ID  object  f(x)  supported  fully  on  the  interval  [—L,L\  has  a  Fourier-series  expansion, 

-.00  pL 

f(x)  —  —  Y"'  Ftexp(in£x/L),  with  —  /  f(x)  exp{— i2n[£  /  (2L)]x}  dx.  (8) 

~L  t  x  '  J~L 

The  Fourier  transform  (FT)  of  this  expansion  yields  the  full  object  spectrum,  F(u),  in  terms 
of  its  spatial- frequency  samples,  F#  =  F{u  =  £/2L),  as 

OO 

F(u)  =  £  Fi  sinc(2wL  —  £),  (9) 

£=—OG 

much  like  the  usual  Shannon-Nyquist  expansion  for  a  band-limited  function  in  terms  of  its 
critical  spatial  samples.  A  similar  but  more  involved  expansion  for  2D  supported  object 
brightness  distribution  may  also  be  given  at  least  for  the  case  of  rectangular  and  circular 
supports,  the  latter  involving  the  Fourier-Bessel  functions  in  place  of  the  Fourier  exponentials 
of  Eq.  (8). 

For  a  one-dimensional  (ID)  imager  with  an  idealized  sensor  array  with  no  dead  space 
between  successive  pixels  and  unit  quantum  efficiency,  the  pixel  transfer  function  (PTF), 
P(u),  has  the  form,  P[u )  =  5w  sine (u5w),  where  the  sine  function  is  defined  as  sinc(a;)  = 
sin  ttx/(ttx).  For  a  clear  aperture,  the  incoherent  OTF  is  (1  —  \u\/B0),  B0  being  the  optical 
bandwidth,  the  blurred  image  data  at  the  kth  pixel  may  be  cast  as  an  integral  over  the 
normalized  frequency,  u  =  u/Bq,  with  all  spatial  variables  scaled  by  the  critical  sampling 
interval,  1/(2F>0), 


X  I 

9k  =  2  J  dU  exP [inu{wk  -t)}  sinc(hy/2)  (1  -  \u\)  F(B0u)  +  nk,  (10) 

where  y  =  2 B0Sw  is  the  ratio  of  detector-sampling  and  critical  optical-sampling  intervals, 
nk  is  any  additive  sensor  noise  at  the  kth  pixel,  and  t  is  any  sub-pixel  shift  for  an  image. 
For  y  <  1,  the  detector  oversamples  the  image  while  for  y  >  1  it  undersamples  the  image. 

When  expansion  (9)  is  substituted  into  Eq.  (10),  one  sees  that  the  image  data  are  sen¬ 
sitive  to  object  Fourier  amplitudes  F(  that  lie  outside  the  optical  bandwidth.  This  is  a 
result  of  support-induced  correlation  of  spatial  frequencies  inside  to  those  outside  the  op¬ 
tical  bandwidth  of  the  imager.  The  Hermitian  FI  matrix  for  the  sensitivity  of  the  image 
data  to  the  spectral  amplitudes  can  be  computed  via  a  quasi- Gaussian  approximation  for 
the  interesting  case  of  Poisson  image  data  in  the  presence  of  Gaussian  additive  sensor  noise. 
The  diagonal  elements  of  the  inverse  of  the  FI  matrix  represent  the  minimum  variance,  also 
known  as  the  Cramer- Rao  lower  bound  (CRB),  with  which  the  real  and  imaginary  parts  of 
the  corresponding  spectral  amplitude  of  the  object  may  be  estimated. 

Figures  8  display  the  results  of  the  CRB  calculation  for  a  ID  Gaussian-shaped  object 
distribution  for  different  data  SNR  values  and  for  the  space-bandwidth  product  (SBP)  Q  = 
2 BqL  =  20.  Note  the  sharp  rise  of  the  minimum  variance  of  estimation  when  spectral  samples 
increasingly  beyond  the  band  edge  are  considered.  Indeed,  the  maximum  extension  of  the 
SBP  that  can  be  achieved  by  means  of  support  knowledge  alone  is  of  order  1,  and  depends 
logarithmically  on  the  SNR. 
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(b)  w  =  0.5,  Peak  SNR  =  10 


Figure  8:  Normalized  CRB  for  a  Fourier  sample  plotted  vs  the  index  of  the  sample  for  Q  —  20 
and  for  two  peak  SNR  values,  103  and  10. 


The  analysis  based  on  the  sampling  expansions  of  the  type  (9)  provides  a  unified  treat¬ 
ment  of  both  digital  (DSR)  and  optical  superresolution  (OSR).  Indeed,  by  analyzing  a  se¬ 
quence  of  sub-pixel- shifted  undersampled  images  one  can  show  the  minimum  variance  needed 
both  to  extract  the  critically  sampled  object  brightness  distribution  and  to  even  achieve 
optical-bandwidth  extension  to  recover  super-critically  sampled,  or  super  resolved,  image 
data.  These  and  2D  image  DSR  and  OSR  considerations  were  presented  in  great  detail  in 
three  papers  (Prasad  and  Luo,  2009  -  all  3  papers  listed  in  Sec.  3). 

2.7  CASSI  Instrumentation 

We  developed  hardware  to  supply  high  resolution,  broad  bandwidth  HSI  data.  This  in¬ 
cluded  a  new  instrument  that  is  capable  of  acquiring  multiple  frames  of  uniquely  coded  data 
over  320nm-700nm.  In  addition,  two  detectors  are  available  for  data  acquisition,  one  for 
UV-visible  and  another  for  higher  frame-rate  visible  data.  Lastly,  a  pushbroom  variant  of 
the  CASSI  system  was  developed  for  ground-truth  data  comparison  with  the  compressively 
sensed  CASSI  data  reconstructions. 

The  new  CASSI  instrument  is  shown  in  figure  9(a).  The  UV-visible  detector  shown  can 
be  swapped  for  a  visible  detector  in  a  matter  of  minutes  depending  on  the  data  acquisition 
needs.  The  UV  detector  has  a  quantum  efficiency  (QE)  of  10%  or  less,  while  the  visible 
detector  has  a  QE  of  55%,  decreasing  integration  time  by  5.5  times.  This  is  particularly 
useful  for  dynamic  scenes  to  decrease  blur  and  increase  frame  rates  for  video  applications, 
snapshots  can  be  grouped  for  much  higher  resolution  reconstructions.  The  relay  lens  was 
designed  by  Duke  and  custom  built  by  Shanghai  Optics. 

The  spectral  range  of  the  new  instrument  is  320-700nm  divided  over  nearly  100  channels. 
The  optics  have  also  been  completely  custom  designed  for  this  specific  application,  greatly 
enhancing  the  performance  of  the  instrument  compared  to  the  previous  generation  system. 
The  camera  is  an  Imperx  Lynx  IPX-4M15  UV-enhanced  2048x2048  monochrome  detector 
(Visible:  Allied  Vision  Technologies,  Pike,  2048x2048).  The  aperture  code  is  a  random 
pattern  lithographically  etched  on  a  chrome  substrate  by  Photo  Sciences  with  minimum 
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(b) 

Figure  9:  (a)  Coded  aperture  spectral  imager;  (b)  custom  relay  lens  design. 


features  corresponding  to  two  pixels  or  14.8 pm.  The  aperture  code  is  imaged  onto  the 
detector  by  the  relay  lens  (custom  fabrication  by  Shanghai  Optics)  and  dispersed  by  a  double 
amici  prism.  The  aperture  code  is  modulated  by  a  piezo  system  (Newport  Corporation)  by 
up  to  21  pixels  on  the  detector.  The  objective  lens  was  provided  by  Coastal  Optics  and  is  a 
true  apochromat  across  the  entire  spectral  range  of  the  instrument. 

At  the  heart  of  the  system  is  a  newly  designed  relay  lens  in  figure  9(b)  consisting  of  a 
modified  Cooke  triplet  with  field  lenses.  The  system  is  symmetric  around  the  prism,  which  is 
in  the  collimated  space  of  the  Cooke  triplet.  The  7  A  pm  pixels  of  the  detector  require  an  MTF 
of  at  least  67.6 £p/mm,  while  the  aperture  code  size  specifies  the  FOV.  The  entire  spectral 
range  of  300  —  700nm  must  be  focused  at  all  wavelengths,  ideally  being  apochromatic.  The 
difficulty  with  this  wavelength  range  is  that  most  glass  types  do  not  transmit  below  350nm. 
Fused  silica  (FS),  calcium  fluoride  ( CaF2 ),  and  magnesium  fluoride  ( MgF2 )  transmit  well 
below  300mn,  but  because  of  MgF2  being  birefringent,  it  is  not  preferred  by  Shanghai  Optics. 
With  only  two  glass  types,  color  correction  is  nearly  impossible,  so  BK7  was  considered, 
which  for  5 mm  thickness  transmits  close  to  90%  above  320 nm.  Compromising  slightly  at 
the  lower  UV  significantly  improved  the  MTF,  with  a  minimum  of  40  —  80 (tp/mm  and  most 
areas  and  wavelengths  above  90 tp/mm. 

Actual  MTF  measurements  of  the  lens  were  done  using  an  USAF  resolution  target,  illu¬ 
minated  by  630nm  (lOnm  bandpass)  light.  An  Aptina  MT9F002  monochrome  14  megapixel 
detector  with  1.4/im  pixels  was  used.  The  Air  Force  target  measures  the  contrast  transfer 
function  (CTF),  which  can  be  converted  to  MTF  approximately  by  multiplying  the  CTF 
by  0.785.  This  is  only  valid  for  higher  frequencies,  shown  in  figure  10(b),  where  frequencies 
below  60 ftp /mm  are  estimated  and  shown  as  the  dashed  line.  This  plot  is  only  for  the  center 
of  the  lens.  The  lOnm  bandpass  region  disperses  about  6.6 pm,  or  4.7  pixels  on  the  Aptina 
detector  in  the  horizontal  direction,  therefore,  only  MTF  measurements  were  made  for  the 
horizontal  bars.  20%  contrast  is  observed  at  lAAip/mm  and  67.6 tp/mm  is  above  50%.  This 
is  in  good  agreement  with  the  MTF  shown  from  Zemax.  Note  in  Zernax  that  the  center 
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Figure  10:  100%  zoomed,  cropped  image  of  1951  USAF  target  illuminated  with  630nm  light 
(lOnm  bandpass)  imaged  on  an  Aptina  MT9F002  monochrome  14MP  detector  with  1.4/zm 
pixels  near  the  center  of  the  field.  There  is  approximately  5  pixels  of  dispersion  in  the 
horizontal  direction;  hence  the  vertical  bars  are  more  blurred. 


MTF  curve  is  worse  than  the  five  and  seven  millimeter  curves  since  it  was  corrected  over  the 
entire  field. 

2.8  Experimental  Results 

Segmentation  in  collaboration  with  WFU  researchers  CASSI  images  of  a  color 
chart  illuminated  using  Solux  daylight  4700K  lamps,  shown  in  figure  11(a)  are  used  for 
reconstructing  both  the  color  segments  and  their  spectra.  Since  the  spectral  signatures 
within  each  segment  are  fairly  uniform,  the  scene  serves  as  a  good  candidate  for  testing 
the  algorithm.  For  brevity,  we  only  show  the  reconstructed  colors  in  the  second  row  as 
indicated  by  the  red  box  in  figure  11(b).  The  spectra  of  five  colors  excluding  the  purple  are 
measured  by  an  Ocean  Optics  (Dunedin,  Florida)  spectrometer.  A  total  of  12  frames  of  SD- 
CASSI  images  are  used  for  reconstructing  a  hyperspectral  cube  in  44  wavelength  channels. 
Figure  11(c)  shows  the  identified  six  segments.  Clearly  we  are  able  to  identify  the  sharp 
boundaries  of  the  segments.  For  the  accuracy  of  reconstructed  spectra,  we  compared  them 
with  both  references  by  the  spectrometer  in  figure  12(a)  and  the  reconstruction  by  TwIST 
in  figure  12(b).  Our  results  match  almost  identically  with  those  by  TwIST  and  show  good 
agreement  with  the  spectrometer  references.  Note  that  there  has  not  been  any  calibration 
to  accommodate  for  quantum  efficiency,  wavelength  offset,  and  channel  non-linearity.  This 
accounts  for  the  discrepancies  shown  in  figure  12(b). 

TwIST  Reconstructions  by  the  Duke  group  The  UV-CASI  system  was  then  taken 
outside  during  an  early  summer  afternoon  when  the  UV-index  was  fairly  high.  Many  flowers 
have  ultraviolet  and  visible  features,  making  an  interesting  target  for  spectral  imaging.  An 
RGB  image  of  an  Ox-Eye  daisy  (Hcliopsis  scabra)  is  shown  below  in  figure  13(a),  with  a 
composite  CASI  image  shown  to  the  right.  The  visible  RGB  image  shows  no  variation  of 
the  petals,  but  solid  yellow  color  throughout.  The  CASI  image  in  figure  13(b)  shows  that 
in  the  UV  bands,  a  dark  region  surrounds  the  central  part  of  the  flower.  Selected  spectral 
channels  are  shown  in  figure  14  from  360nm  to  702nm. 
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Figure  11:  (a)  RGB  image  of  the  Gretagmacbeth  Color  Checker,  (b)  The  original  SD-CASSI 
image  with  the  reconstructed  area  in  the  red  rectangular  box.  (c)  The  reconstructed  seg¬ 
ments. 


Two  separate  images  were  acquired  with  the  CASI  system  since  there  are  no  bandpass 
filters  available  in  the  range  of  320-700nm  that  strongly  filter  the  IR  bands.  The  visible  band 
was  imaged  using  a  Baader  UV/IR-cut  filter,  blocking  below  400nm  and  above  680nm.  The 
UV  bands  were  imaged  using  a  Badder  U-hlter  Venus,  blocking  both  the  visible  and  IR.  The 
Imperx  camera  is  sensitive  between  200nm  and  lOOOnm,  so  blocking  the  700-1000nm  region 
is  necessary.  The  reconstructions  shown  in  figure  14  were  done  using  TwIST  and  20  uniquely 
coded  frames  for  UV,  and  24  frames  for  the  visible. 

Spectral  accuracy  Spectral  features  for  a  non-spatially  varying  scene  are  shown  in  figure 
15.  An  Ocean  Optics  USB2000  spectrometer  is  compared  to  the  CASSI  system  with  a  1- 
pixel,  slit  and  a  random  mask  with  2-pixel  features.  The  slit  and  random  mask  have  similar 
spectra,  with  slight  offsets  from  the  ocean  optics.  The  shown  spectra  are  uncalibrated  for  the 
non-uniform  quantum  efficiency  of  the  detector  and  optical  transmission.  The  slit  spectra 
for  CASSI  is  directly  read  from  the  detector  cross-section,  while  the  random  code  data  was 
reconstructed  using  TwIST  as  in  figure  13(b).  The  random  code  has  slightly  better  resolution 
than  the  slit,  even  though  the  slit  is  narrower  than  the  minimum  feature  of  the  random  code. 


Pushbroom  versus  CASSI  The  coded  aperture  system  can  be  compared  to  a  reference 
data  cube  by  swapping  the  coded  aperture  with  a  vertical  slit  and  pushbrooming  across  the 
detector.  The  results  for  this  experiment  are  shown  in  figure  16. 
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Figure  12:  (a)  The  reconstructed  spectra  compared  with  reference  by  an  Ocean  Optics 
spectrometer,  (b)  The  reconstructed  spectra  compared  with  reconstruction  using  TwIST. 


Figure  13:  (a)  Visible  image  taken  with  standard  RGB  camera,  (b)  CASSI  image  with  UV 
channels  on  the  left  and  visible  on  the  right. 


2.9  Joint  segmentation  and  reconstruction  methods 

We  developed  numerical  methods  for  the  joint  reconstruction  and  segmentation  of  spectral 
images  taken  by  compressive  sensing  coded  aperture  snapshot  spectral  imagers  (CASSI) 
(Zhang,  Plemmons,  Kittle,  Brady  and  Prasad,  2011).  In  a  snapshot,  a  CASSI  captures 
a  two-dimensional  (2D)  array  of  measurements  that  is  an  encoded  representation  of  both 
spectral  information  and  2D  spatial  information  of  a  scene,  resulting  in  significant  savings 
in  acquisition  time  and  data  storage.  The  reconstruction  process  decodes  the  2D  measure¬ 
ments  to  render  a  three-dimensional  spatiospectral  estimate  of  the  scene,  and  is  therefore  an 
indispensable  component  of  the  spectral  imager.  Because  often  the  scene  to  be  reconstructed 
from  the  HSI  data  typically  consists  of  spectrally  and  spatially  homogeneous  segments  that 
can  be  represented  sparsely  in  an  appropriate  basis,  as  is  generally  the  case  with  SOI,  we 
seek  a  particular  form  of  the  compressed  sensing  solution  that  assumes  spectrally  homoge¬ 
neous  segments  in  the  two  spatial  dimensions,  and  greatly  reduces  the  number  of  unknowns, 
often  turning  the  under- determined  reconstruction  problem  into  one  that  is  over-determined. 
Numerical  tests  are  reported  on  both  simulated  and  real  data  representing  compressed  mea¬ 
surements. 
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As  an  illustration  we  solve  the  joint  reconstruction/segmentation  problem  by  integrating 
a  spectral  signature  solver  with  a  segmentation  solver  -  our  fuzzy  segmentation  method 
for  hyperspectral  images.  A  sample  result  from  this  study,  using  a  simulated  hyperspectral 
image  of  the  Hubble  Space  Telescope  satellite  developed  in  our  papers  (Gillis  and  Plemmons, 
2011;  Zhang,  Wang,  Plemmons,  and  Pauca,  2008),  is  given  in  Figure  17.  Eight  materials 
typical  to  those  associated  with  satellites  are  used,  Hubble  aluminum,  Hubble  glue,  Hubble 
honeycomb  top,  Hubble  honeycomb  side,  solar  cell,  bolts,  rubber  edge,  and  copper  stripping. 
Note  the  excellent  recovery  of  six  of  the  seven  material  spectra;  the  discrepancy  between  the 
estimated  and  true  spectra  for  the  bolts  segment  is  related  to  its  relatively  small  abundance 
in  the  object. 

2.10  Coupled  segmentation  and  denoising/deblurring  methods 

A  crucial  aspect  of  spectral  image  analysis  is  the  identification  of  the  materials  present  in 
the  object  or  scene  being  imaged  and  to  quantify  their  fractional  abundances  in  the  mixture. 
An  increasingly  useful  approach  to  extracting  such  underlying  structure  is  to  employ  image 
classification  and  object  identification  techniques  to  compressively  represent  the  original  data 
cubes  by  a  set  of  spatially  orthogonal  bases  and  a  set  of  spectral  signatures.  Owing  to  the 
increasing  quantity  of  data  usually  encountered  in  hyperspectral  data  sets,  effective  data 
compressive  representation  is  an  important  consideration,  and  noise  and  blur  can  present 
data  analysis  problems.  We  developed  image  segmentation  methods  for  hyperspectral  space 
object  material  identification  and  coupled  the  segmentation  with  a  hyperspectral  image  data 
denoising/deblurring  model  to  propose  this  method  as  an  alternative  to  a  tensor  factorization 
methods  proposed  recently  (Zhang,  Wang,  Plemmons,  and  Pauca,  2008).  The  model  provides 
the  segmentation  result  and  the  restored  image  simultaneously.  Numerical  results  show  the 
effectiveness  of  our  proposed  combined  model  in  hyperspectral  material  identification  (Li, 
Ng,  and  Plemmons,  2011). 

2.11  Sparse  Nonnegative  Underapproximation  for  Hyperspectral 
Image  Analysis 

Dimensionality  reduction  techniques  such  as  principal  component  analysis  (PCA)  are  power¬ 
ful  tools  for  the  analysis  of  high-dimensional  data.  In  hyperspectral  image  analysis,  nonneg¬ 
ativity  of  the  data  can  be  taken  into  account,  leading  to  an  additive  linear  model  called  non¬ 
negative  matrix  factorization  (NMF),  which  improves  interpretability  of  the  decomposition. 
Recently,  another  technique  based  on  underapproximations  (NMU)  has  been  introduced, 
which  allows  the  extraction  of  features  in  a  recursive  way,  such  as  PCA,  but  preserving 
nonnegativity,  such  as  NMF.  Moreover,  in  some  situations,  NMU  is  able  to  detect  automat¬ 
ically  the  materials  present  in  the  scene  being  imaged.  However,  for  difficult  hyperspectral 
datasets,  NMU  can  mix  some  materials  together,  and  is  therefore  not  able  to  separate  all  of 
them  properly.  In  this  paper  we  introduce  sparse  NMLT  by  adding  a  sparsity  constraint  on 
the  abundance  matrix  and  use  it  to  extract  materials  individually  in  a  more  efficient  way 
than  NMU.  This  is  experimentally  demonstrated  on  the  HYDICE  images  of  the  San  Diego 
airport  and  the  Urban  dataset  (Gillis  and  Plemmons,  2011). 
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Figure  14:  UV-Visible  reconstruction  from  360nm  to  700nm,  taken  using  multiple  frames  in 
direct  sunlight. 
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Figure  15:  Comparison  of  the  CASSI  spectral  resolution  with  a  1-pixcl  wide  slit,  2-pixel  fea¬ 
ture  random  mask  reconstructed  with  TwIST  and  24-frames,  and  commercial  Ocean  Optics 
spectrometer.  Both  CASSI  spectra  are  uncalibrated  for  quantum  efficiency  of  the  detector, 
optical  transmission  efficiency,  and  channel  sensitivity. 
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Figure  16:  Comparison  with  a  commercial  Nikon  SLR  camera  in  monochrome,  CASSI  with 
the  coded  aperture  replaced  by  a  slit  and  pushbroomed,  and  a  multiframe  CASSI  image. 
Acquisition  time  for  the  slit  was  around  20  minutes,  for  CASSI  about  5  seconds. 
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Figure  17:  (a)  Raw  CASSI  simulated  compressed  sensing  image.  (b)  The  recon¬ 

structed/segmented  CASSI  image.  False  color  identifies  the  materials  in  the  satellite,  (c) 
The  estimated  material  spectral  signatures  compared  with  the  original  ones. 
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